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We study the dynamics of elastic interfaces — membranes — immersed in thermally excited fluids. 
The work contains three components: the development of a numerical method, a purely theoret- 
ical approach, and numerical simulation. In developing a numerical method, we first discuss the 
dynamical coupling between the interface and the surrounding fluids. An argument is then pre- 
sented that generalizes the single-relaxation time lattice-Boltzmann method for the simulation of 
hydrodynamic interfaces to include the elastic properties of the boundary. The implementation of 
the new method is outlined and it is tested by simulating the static behavior of spherical bubbles 
and the dynamics of bending waves. By means of the fluctuation-dissipation theorem we recover 
analytically the equilibrium frequency power spectrum of thermally fluctuating membranes and the 
correlation function of the excitations. Also, the non-equilibrium scaling properties of the mem- 
brane roughening are deduced, leading us to formulate a scaling law describing the interface growth, 
W 2 (L,T) = L 3 g(t/L 5 / 2 ), where W, L and T are the width of the interface, the linear size of the 
system and the temperature respectively, and g is a scaling function. Finally, the phenomenology 
of thermally fluctuating membranes is simulated and the frequency power spectrum is recovered, 
confirming the decay of the correlation function of the fluctuations. As a further numerical study 
of fluctuating elastic interfaces, the non-equilibrium regime is reproduced by initializing the system 
as an interface immersed in thermally pre-excited fluids. 



I. INTRODUCTION 

Elastic interfaces in fluids, such as biological mem- 
branes, have spurred a strong interest in recent years 
not only because of their practical importance, but also 
because of the intriguing complexity of their phenomenol- 
ogy. Many aspects of this rich subject have been stud- 
ied from different approaches. Since temperature and its 
effects play a primary role in natural phenomena, ther- 
mal fluctuations of fluid, hexatic, nematic and polymer- 
ized membranes have been theoretically discussed in a 
number of works |l]-|(|. Also, due to their relevance to 
mechanical and chemical interactions in biological sys- 
tems, shape transformations and fluctuating topologies of 
elastic interfaces immersed in fluids with, in many cases, 
shear flows have been thoroughly studied [0-^2). Crys- 
talline membranes, also known as polymerized or teth- 
ered membranes, are, in particular, a fascinating subject 
with important concrete realizations in nature, such as 
the cytoskcleton of mammalian erythrocytes (red blood 
cells) Jl^,|l4| • Among the systems that can be physically 
realized in a laboratory, inorganic crystalline membranes 
were examined in [[15 16 1. On the other hand, theoreti- 
cal work on tethered membranes has proved successful in 
addressing the crumpling transition jT7|-pl|. Other the- 
oretical results on finite-size effects in fluid membranes 
can be found in |22U24|. Furthermore, some books and 



reports p5|-p9[ reviewing the statistical mechanics, ther- 
modynamics and geometrical structure of membranes 
have shed further light on the understanding of the sub- 
ject. 

In the present work, we focus on the physics of elastic 
interfaces in fluids when thermal fluctuations in the bulk 
generate correlated forces and subsequent excitations on 
the boundary. This dynamical coupling between inter- 
face and surrounding fluid is believed to be responsible 
for interesting phenomena, such as the flickering of ery- 
throcytes Q or the physical distribution of particulates 
inside certain lipid vescicles [^0|. Because many such 
problems do not have a closed-form solution, numeri- 
cal simulations can provide us with a valuable tool for a 
deeper understanding and general guidelines for design- 
ing new experiments. More generally, the ability to sim- 
ulate thermally fluctuating elastic membranes in fluids 
allows a computational model to be effectively employed 
in the simulation of the biophysical systems for which 
the effects of a finite temperature are to be accounted 
for. Computational works on fluctuating membranes 
have appeared rarely in literature [fn 31 - 3j| . As an inter- 
esting example, in the work by Goctz et al. fl3l| the power 
spectrum of a fluctuating bilayer membrane in vacuum 
is obtained by molecular dynamics. Also, the reader 
may find in |3^] an example of how a mean-field theory 
approach in the framework of a lattice-Boltzmann model 
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can be effectively employed in the study of phase sep- 
aration with boundaries driven by surface tension and 
bending stiffness. Broadly speaking, elastic forces are 
governed by the local curvature of the interface and their 
wavelengths are two to four orders of magnitude larger 
than molecular ones. Therefore, a major difficulty in 
modeling such systems is one of describing the problem 
at the different length and time scales of molecular and 
elastic interactions in an unified and consistent approach. 

Our work is threefold. First, we create a new numeri- 
cal method for the study of fluctuating membranes. Sec- 
ond, we use the method to simulate phenomena associ- 
ated with the coupling of the interface and the surround- 
ing fluids. Third, after reviewing the properties of fluc- 
tuating membranes at equilibrium, we deduce the non- 
equilibrium scaling law governing the interface roughen- 
ing and show that our predictions are observed in simu- 
lations. 

We develop a lattice-Boltzmann model. Not only does 
our method allow the simulation of bending waves and 
fluctuating membranes in fluids, but it also allows the 
study of more complicated physical problems, in which 
the interface has many distinct components and the fluids 
have prescribed velocity vector fields, inducing stresses 
on the boundary. So that our method may eventually be 
used to simulate complex flows in complex geometries, it 
is designed to produce thin interfaces, with a thickness 
of the order of few (~ 3) lattice units. Accordingly, we 
choose to not describe the elastic boundaries by means of 
a slowly varying order parameter. We choose instead to 
represent the dynamics of interfaces with bending rigid- 
ity by means of a free energy |l| |I"l] , |3~l ■ in which the loca- 
tion and the geometrical properties — the curvature — of 
the membrane appear explicitly. Thermal fluctuations 
are introduced in the model as a Gaussian noise in the 
lattice-Boltzmann equation 



35 36 1 . The link between dif- 
ferential geometry and microscopic dynamics is provided 
by a suitable perturbation, driven by the local curva- 
ture, of the occupation numbers, similarly to a procedure 
already successful in the study of interfaces with surface 
tension j37|. 

This paper is organized as follows. In Section || we 
present an overview of the relevant interface and fluid 
dynamics. Section III focuses on the fluid-interface cou- 



and compare them with the theo ry previously outlined. 
Conclusions follow in Section VII. 



pling and on how the macroscopic equations of motion for 
the membrane are translated into microscopic mechanical 
prescriptions for the occupation numbers. We thus pro- 
vide the theoretical basis for building a lattice-Boltzmann 
computational model. Section IV reports the results of 
the simulation of spherical bubbles and bending waves. 
Also, the experimental dispersion relation is here com- 
pared to the theoretical prediction. In the following sec- 
tion we then discuss the theoretical description of the 
physics of one-dimensional fluctuating membranes cou- 
pled to thermally excited fluids. In Section [VI 



we use our 



model to simulate fluctuating elastic interfaces and study 
both the non-equilibrium roughening and the stationary 
state. We present here the results of our computations 



II. DYNAMICAL COUPLING OF ELASTIC 
INTERFACES TO AN EXTERNAL FORCE 

In this section we provide a theoretical motivation for 
the lattice-Boltzmann microscopic dynamics that con- 
stitutes the basis of our model. We recall the com- 
mon expression of the free energy for elastic interfaces 
and apply Hamilton's variational principle to recover the 
macroscopic equation of motion. 

The dynamics of membranes with bending or flexural 
rigidity e and surface tension a is assumed to be governed, 
for the case of a vanishing spontaneous curvature, by the 
free energy [EJfll 34 



T 
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e / dSH 2 + e / dSK + a dS, 



(2.1) 



where the integrations are performed over the area of a 
2D membrane or over the length of a ID interface. H is 
the sum of the principal curvatures and K is their prod- 
uct (the Gaussian curvature) . e is named the saddle-splay 
modulus and the second term in the above equation is 
responsible for the energy gain/loss due to a change in the 
interface genus, according to the Gaufi-Bonnet theorem 
J dSK — 2irx = &ir(l — g), where \ is the Euler-Poincare 
characteristic of the surface and g is its genus pg[ ]. For 
the ID interfaces studied here, this term acts as a spon- 
taneous curvature. Since we set it to vanish, the saddle- 
splay term will be neglected. The dynamical coupling of 
the membrane with its surroundings is realized by intro- 
ducing the force F per unit area or unit length exerted 
by the fluids on the interface. By appl yin g Hamilton's 
variational principle to the free energy (pj) and includ- 
ing a term for the work done by F when the membrane 
undergoes a configurational change [p] j39| , ^0| | , the follow- 
ing relation between F± , the component of the force per- 
pendicular to the interface, and the geometric properties 
is derived in Appendix |X] 

F± = oH + eHj2"/f -eV 2 i? . (2.2) 

Here the 7^'s are the principal curvatures of the (n — 1)- 
dimensional hypersurface. This expression specializes to 
the case of a ID membrane as 
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(2.3) 



In ( |2.3| ) s is the arc length of a canonical parameterization 
of the interface and 7(5) is the local curvature. 



III. NUMERICAL METHOD 



The lattice-Boltzmann method [ p7|j4l| ] we adopt solves 
the incompressible Navier-Stokes equations for the fluid 
dynamics, namely 
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Here u is the fluid macroscopic velocity, p is the local 
pressure, p the density and p, is the viscosity coefficient. 
We show in this section how u, p and p are expressed 
in terms of microscopic quantities and how the cou- 
pling betw een the fluid dynamics (3.1) and the interface 
dynamics (2.3) is realized by means of the microscopic 
pressure tensor. 

In making the connection between the macroscopic 
relation (2.3) and the microscopic dynamics of our 
lattice-Boltzmann method, we shall follow a procedure 
that has been previously employed for interfaces with 
surface tension |}7|,[l2). Notice first that F± corresponds 
to the local fluid pressure gap Ap = p\ — pi across the 
interface, so that one can rewrite (|2.3|) as 



(3.2) 



Ap(s) — cr'-f(s) + e-f 3 (s) — £7(s), 



where a dot denotes the derivative with respect to the 
arclength. The mechanical relation |57 43 between the 
normal pressure and the longitudinal pressure p t reads 



Ap(s)=j(s) \p„(s) -p t (s,y)]dy 



(3.3) 



where the integration is carried out along the y direction 
perpendicular to the interface. The average between the 
pressures at either side of the membrane is p n = (pi + 
p 2 )/2 = P2 + Ap/2. In practice Ap < (p n -p t ) < p n . We 
therefore si mply replace p n (s) withp„(s) «p 1 (s) ~ P2{s) 
and recast (3.2) as 



7 / [Pn- Pt] dy = ct7 + £7 3 



£7. 



(3.4) 



In the lattice-Boltzmann method j3^ , |4l| ], the momen- 
tum flux tensor II is expressed in terms of the discretized 
velocity distribution functions m(x.,t) 



n(x, t) = m(x, t)ciCi . 



(3.5) 



where rii(x, t) represents the positive real- valued occupa- 
tion number at a given site x and time t with velocity . 
Similar expressions hold for the fluid density and velocity 



P = ^2n l (x,t) 

i 

pu = y^nj(x, t)cj 



(3.6) 
(3.7) 



In the present work, interfaces separate the two compo- 
nents, called "red" and "blue", of a binary fluid. The 
sum of the red occupation numbers fj and the blue ones 
hi at each lattice site is preserved: 



Starting from this distinction in terms of color attributes, 
it is possible to build a numerical method and study mis- 
cible and immiscible fluids, and, after suitable general- 
izations, the physics of even more complex fluids can be 
discussed ^Q. 

The different components of II are related to the 
macroscopic pressure tensor. Let us assume, for exam- 
ple, that the interface is oriented parallel to the x-axis in 
a neighborhood of a given point (x v , 0). Then IL^ and 



n xx replace p n and p t respectively in (3.4), resulting in 



^2 ^2n t (x p ,k,t) 



,)=a 7 + £ 7 3 - £ 7 , (3.9) 
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where k is the discrete coordinate running along the y- 
axis. 

The microscopic dynamics governing the time evolu- 
tion of the distribution functions is described by the 
Boltzmann equation. Lattice-Boltzmann models rely on 
a discretized form of it, namely p7 41 45|] 



n,(x + Ci,t + 1) = ni(x,i) + A 4 [n(x, t)], 



(3.10) 



where the last term is the collision operator, which 
accounts for the change in the occupation numbers due 
to collisions at the lattice sites. For practical reasons, 
Aj[n(x, t)] is usually replaced by its linear expansion 
around the equilibrium populations ^ q (x, t) |ff6|| 

Ai[n(x,t)]=J2 L H h(x,t)-n° q (x,i)], (3.11) 



where Ly is a matrix of costant coefficients. A further 
simplification consists in substituting Lij with a diago- 
nal operator Xs^ij-, where As is the relevant eigen value 
of the linearized Boltzmann operator, so that (|3.10| ) sim- 
plifies in the single relaxation time lattice-Boltzmann 
model 



rii(x + Ci,t + 1) = (1 + X B )ni(x, t) ~ X B n i q (x, t) . 

(3.12) 

This expression conserves mass and momentum. Also, 
the populations rij, in the absence of external forcing, 
converge to the equilibrium occupation numbers n^ q . 

Thermal fluctuations are introduced in the model by 
ad ding a stochastic term A[ (x, t) to the right-hand side 
of fl3.12D such that || 



Aj(x,t) K^^(x,i) (a^ - c 2 S aP /D) . (3.13) 

Here D is the dimensionality of the lattice and the 
random fluctuations a' a g are uncorrelated in space and 
time , and sampled from a Gaussian distribution such 
that 



ni(x, t) = rj(x, t) + bt(x,t) . 



(3.8) 



(a^(x,i)a; c (x',0) = 
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(3.14) 



where the variance A is related to the effective tempera- 
ture T of the fluid. 



A = 2pvk B T\% , 



(3.15) 



by means of the fluctuation-dissipation theorem [ p5[ . 

In order to reproduce the desired surface tension and 
bending stiff ness of the interface, so that the left-hand 
side of (3.9) does not vanish in the proximity of the 
boundary, one may suitably perturb the single relaxation 
time lattice-Boltzmann model by adding a ter m A" to 



the right-hand side. Our perturbation of ( 3.12 ) reads 



A'/ = (S + E 7 2 - £7/7) |f I {dcCip 

a0 



C d a f3/D) -jg- 
(3.16) 



where S and E are two adjustable parameters corre- 
sponding to the physical surface tension a and bending 
rigidity e, and f (x, f) is the local color gradient |3^] as 
defined in Appendix |b[ In order to measure the geo- 
metric properties of the membrane — the curv ature 7 and 
its derivatives — appearing in equation ( 3.16 ), we adopt 
an explicit procedure of localization of the different con- 
nected components of the interface. The method is 
described in some detail in Appendix [§]. There we show 
how, for each component, we map the boundary piece- 
wise to polynomials, from which we extract the curvature 
and its derivatives. This process does not result in a sig- 
nificant time consumption for the whole simulation, as it 
is performed only once at the beginning of each time step 
and it is limited to the lattice sites which constitute the 

interface. 

Expression ( 3.16| ) preserves the total mass and momen- 
tum at a given lattice site, as one can verify by recalling 
the identities 

c ia = 

i 

^ c ia c ll3 = b m c 2 8 a fi/D 

i 

i 

which hold true for tensors with hypercubic symmetry. 
By generalizing the argument set forth in Se ction 10.2 
of (37]], it can be show n th at the inclusion of ( 3.16 ) into 
the right-hand side of (3.12) generates a surface tension a 
and a bending rigidity e of the interface, which are related 
to S and E respectively through the linear relation 
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(3.17) 



valid for the face-centered hypercubic lattice (FCHC), 
which we shall employ throughout this work. The above 



result follows from replacing equations (10.14) and (10.3) 
in |l37ll with o ur eq uati ons (O) ) and the result of the per- 



3.16) and carrying out an analysis 



turbation of ( 3.12| ) by 
similar to that one in 



IV. TESTING THE MODEL 

The surface tension parameter S is set to zero s o that 
purely elastic effects can be studied. In section IV A, 



the results for a spherical bubble surr ounde d by fluids 
on both sides are presented. In section IV B| , the bend- 
ing wave dispersion relation is discussed by studying the 
damped oscillations of sinusoidal interfaces. 



A. Spherical bubbles 

The lattice-Boltzmann simulation is initialized as a 
spherical bubble of a fluid, here called "red" for prac- 
tical purposes, immersed in a bath of "blue" fluid with 
linear dimension L much larger than the radius R of the 
sphere. The pressure gap Ap between the red inner part 
of the bu bble and the outside blue sea is predicted by 
replacing (3.17) in (3.2), resulting in 



Ap = 



192p£_ 
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(4.1) 



as a is set to vanish in our experiments, 7=1 / R and 7 = 
for a sphere. The pressure is computed by measuring 
the fluid density according to the equation of state (in 
the absence of a net momentum) [§7| 



P(P) 



b m p 
2b ' 



(4.2) 



where b m is the number of velocity vectors (24 for 
FCHC), while b = b m + b r includes the number b r of 
rest particles (in this case b r = 16). Simulations were 
performed for different bubble sizes, ranging from R = 8 
to R = 64, with average density of 0.5 particles per lat- 
tice site, Xb = — 1 and E = 10~ 3 . The measured values 
of Ap are plotted versus R~ 3 in Figure 0. The agreement 
with the predicted relation Ap = 0.096i?~ 3 , drawn as a 
solid line, is very good. For bubble radii smaller than 
eight lattice units, discretization effects cause the mea- 
sured curvature to be off more than 15%, and ultimately 
wrong when the radius of curvature is of the order of the 
interface thickness, that is about four lattice units. 



B. Bending waves 

While bubble pressure gaps verify the equilibrium 
properties of the membrane, the simulation of bending 
waves provide an effective tool for testing the dynamics 
of the fluid-interface coupling. In testing the bending 



4 



x 10 




FIG. 1. Verification of the Ap 



law for a fluid bubble 



immersed in an immiscible sea. The solid line is the theoret- 
ical prediction for the parameters specified in the text. The 
pressure gap Ap is determined from the equation of state 
p(p) = 3p/10, where p is the fluid density. The radius R of 
the bubble is expressed in lattice units. 
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FIG. 2. We report in figure the results of the lat- 
tice-Boltzmann simulation of a bending wave with wavelength 
A = 100 lattice units (circles). The theoret ical predic tion s 
from the normal- mode analysis, equations (O) and (4.4), 
are graphed as a solid line. hk 1 (t) is the Fourier transform 
in lattice units of the interface profile, corresponding to the 
wavenumber ki — 2n/X. 



wave dispersion relation, we initialize the system as a 
square region filled with two immiscible fluids with the 
same densities and viscosities. Such fluids are separated 
by a sinusoidal interface, whose wavelength A is equal to 
the linear dimension of the box. We impose periodicity 
along the horizontal axis, while free-slip boundary condi- 
tions are prescribed along the vertical axis. The damped 
oscillations of bending waves for membranes immersed 
in viscous fluids have been analytically studied in the lit- 
erature [j50| . For the initial conditions we impose one 
expects the time-dependence of the first normal mode to 
be described by 



h kl (t) - h° kl cos[3t(u;)t] e 



(4.3) 



where k\ = 2n/L, with L the size of the box, and h® is 
the initial amplitude of the sinusoidal wave. The complex 
angular frequency u> is related to the wave number k = 
2tt/X through the dispersion relation H 



o 2 = — (l-k/q) 
q = \/k 2 — iu) jv 



(4.4) 



where v = fi/p is the kynematic shear viscosity. The 
numerical simulation of a bending wave for a sys- 
tem defined by E = 0.25, L = 100, p = 0.5, \ B = 
— 1.0, h? ki / L — 0.05 is giv en as an example in Figure ||. 
By means of (3.17), (4.4) and the relation between the 
viscosity v and the Boltzmann eigenvalue Xb [|37|, 
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(4.5) 



Therefore one may fully predict the behavior of the nor- 
mal mode ( |4.3| ) . hj Cl (t) was computed as the Fourier 
transform of the experimental interface profile h(x,t), 
recorded at each time step . Figure || shows the evolution 
of hk 1 (t) according to ( [4.3| ) (solid line) and our numerical 
results (circles). By fitting the first c ycle of the time evo- 
lution of /ifcj (t) to a curve of the form (4.3) is initially 
prescribed) we collected the numerical data about the 
complex angular frequency uj, whose real and imaginary 
parts correspond to the oscillation frequency and damp- 
ing rate, respectively, of the bending wave. We repeated 
the simulation of Figure ||for different wavelengths, col- 
lecting the data about the damping rates 3(w) and the 
oscillation frequencies In Figure]^ we produce the 

experim enta l results together with the numerical solu- 
tion of (4.4), We notice that the agreement between 
experimental data and theoretical predictions improves 
for larger systems and the relative errors reduce to few 
percent for wavelengths of 0(1O 2 ) lattice units. The rel- 
ative error is here defined in terms of the theoretical and 
experimental time evolution of the damped wave as 



\ 



hk u ex{t)/hk u th{t)Y 



(4.6) 



t=o 



where r is twice the oscillation period of the bending 
wave. Table [| reports the relative errors for different 
wavelengths, confirming their fast convergence towards 
0(1O -2 ) when A > 70. The error at small wavelengths 
is due to a numerical artifact that manifests itself as an 
effective surface tension. We po stpone further discussion 
of this subject to Section VI B. 
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FIG. 3. (a): Oscillation frequencies collected from the sim- 
ulation of bending waves with wavelengths ranging from 50 
to 140 lattice units are rep resented with filled circles. The 
numerical solution of (1.4) is drawn as a solid line, (b): 
Damping rates for the same experiment correspond to the 
imaginary part of the complex frequency Sf(a>). The discrep- 
ancy at shorter wavelenghts is due to the dampi ng ac tion of 
the effective surface tension discussed in Section VI B. 



TABLE I. Relative erro rs e(A ) for the sinusoidal bending 
waves discussed in Section IV B and in Figure p|. 



A 



£ (A) 



50 
60 
70 
80 
90 
100 
110 
120 
130 
140 



0.2741 
0.1915 
0.1262 
0.0833 
0.0640 
0.0519 
0.0428 
0.0371 
0.0323 
0.0289 



V. THERMAL FLUCTUATIONS OF ELASTIC 
INTERFACES: THEORY 

We now turn to a study of fluctuating membranes. 
This section introduces the theoretical results relevant 
to our discussion. A comparison with numerical experi- 
ments will then be presented in succeeding sections. The 
fluctuation-dissipation theorem provides us with a unified 
description of the steady state and the non-equilibrium 
growth (roughening) of fluctuating interfaces. Emphasis 
is given to the frequency power-spectrum, as it conveys 
all the relevant information about the decay of the cor- 
relation function of the fluctuations. 



A. Correlation functions from the 
fluctuation-dissipation theorem 

A detailed discussion of the fluctuation-dissipation the- 
orem and fluctuating hydrodynamic interfaces driven by 
surface tension has been presented in |5l],[52). Here we 
summarize the theoretical results for ID interfaces with 
bending stiffness. 

The interface height will be denoted by h(x,t) and its 
Fourier transform in wave- number space by h k (t). We 
also introduce the Fourier transform in time as 



h k (u) 



dt h k (t) e 



(5.1) 



where Q is the size of the time integration domain. When 
— > oo the integral in (5.1) will be replaced by ^ — -> 
If 

2-7T J ' 



It can be shown |52|]55| ] that the frequency power- 
spectrum can be written in terms of the response function 
Tfe(w) as 



2ir z uj z L 



where 



2top 



ik(l - k/q) 



efc 4 



(5.2) 



(5.3) 



and q is given in (4.4). The frequency power-spectrum 
carries information about the fluctuation correlation 
functions. In deed, by means of the Wiener-Khintchine 
relation, (5^2) can be used to show that 



(hk(t)hl(0)) 



k B T 
2ttL 



did 



- l ^T k {uo). (5.4) 



The integration in (5.4) can be performed in the complex 
w-plane by means of the contour drawn in Figure |4|, so 
that 



(h k (t)hl(Q)) = (h k (t)h* k (0)) bw + (h k (t)h* k (0)) cut (5.5) 
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FIG. 4. Integration contour for the fluctuation correlation 
function (5.4). The poles in the lower half-plane corresponds 
to the bending wave contribution, while the pole at the origin 
represents the steady state power-spectrum. 



is the sum of the contributions from the two poles (bend- 
ing wave) and from the branch cut. In the case of oscil- 
lating bending waves the contribution from the two poles 
results in 



(h k (t)h* k (0)) 



bw 



k B T 
ek^L 



cos|K(u;)t| 



SR(w) 



sin|9ft(o;)i| 



(5.6) 



where uj as a function of k is given by the dispersion rela- 
tion (4.4). Also, in the long time limit, the contribution 
from the branch cut reads 



<M*K(o)) cut = 



2 k 8 L 



\t\ 3 / 2 



(5.7) 



One can show by means of the fluctuation-dissipation 
theorem 0| that the correlations among the fluctuating 
forces acting on the interface, and due to the uncorrelated 
thermal excitations in the bulk, have indeed the same 
te mpor al decay as ( |5.7P . A special case is the one of t = 
in (5.4). The integration contour can then be closed in 
the upper half-plane, including only the pole at uj = 0. 
As a result 



h k (0)\' 



k B T 
ek*L 



(5.8) 



represents the mean square amplitude of the fc-th 
mode [|| , which can also be understood by applying the 
energy equipartition theorem to the free energy (2.1). 



B. Interface roughening 

The mean-square width W 2 (L,t) of an interface with 
vanishing mean height, defined as 



W 2 (L,t) = \ [ dxh 2 (x,t) 
L Jo 



(5.9) 



will here be used to describe th e in terface roughen- 
ing. By using Parseval's relation, (5J5) can be recast as 
W 2 (L, t) = J2 k \hk(t)\ . At the steady state, the average 
width can be evaluated by means of (|5.8|), 



W\L) = 



E 



71=1 V 7 



k B T 
720 e 



(5.10) 



where the factor 2 accounts for the two possible orien- 
tations (for a 1-dimensional interface) of the wave vec- 
tors. This result shows that at equilibrium W/L oc L 1 / 2 , 
that is the relative width apparently increases indefinitely 
with the square root of the linear size of the computa- 
tional box. 

In reality, the argument given above and culminating 
with equations ( |5.8| ) and ( |5.10 ) implicitly approximate 
the free energy T — | e J ds j z with the more convenient 
expression T &v — dx (d 2 h / dx 2 ) 2 . Such approxima- 
tion is certainly plausible when 7 w d 2 h/dx 2 , that is 
when the interface-width/lcngth ratio is small. In prac- 
tice, however, when W/L > 1/2 the local curvature is 
significantly different fr om its linearization d 2 h/dx 2 , and 
the original expression (2.1) for the free energy should be 
considered. Therefore (5.8) and ( 5.10D hold true only for 
small values of W/L, while when W/L > 1/2 nonlinear 
terms in the curvature and consequently in the dynam- 
ics prevent the interface-width/lcngth ratio from growing 
indefinitely with the size of the system. On the other 
hand, the lattice-Boltzmann method we are currently 
describing does not suffer from such a limitation. Indeed, 
the interface dynamics discussed in S ectio n || originates 
in our model from the perturbation ( 3.16| ), which corre- 
sponds to the exact analytical form of the elastic force 
F c i = —ed 2 j/ds 2 and not to its linearized approximation 
F c \ ps —ed 4 h/dx i , deducible from !F a p. This allows us to 
study and simulate the dynamics of membranes in the 
nonlinear regime. 

In order to analyze the non-equilibrium roughening of 
the interface, two immiscible fluids at the same tempera- 
ture T are brought in contact at time t = 0. The interface 
between them is initially flat. Due to the thermal fluctu- 
ations in the bulk, standing bending waves will be excited 
on the interface with frequencies given by the dispersion 
relation 



uj [k) 




(5.11) 



wh ere f or simplicity we consider the inviscid limit (y = 0) 
of (4.4). In analogy with the argument of [ pT]j52] ] and by 
means of the fluctuation-dissipation theorem, one antic- 
ipates that the power spectrum is given by 



h k {t)\- 



2k B T 
eLk A 



sin 2 [w (fc)t] 



(5.12) 



so that the mean-square width of the interface is expected 
to grow as 
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W 2 {L,t) 
^ V 27m / 



Yp \T) 



5/2 



(5.13) 



where the thermal average of W 2 i s und erstood. The 
inviscid form of W 2 (L, t) as given by (5.13) does not hold 
for long times, since it does not relax to the equilibrium 
( |5.10 ) as it should if the effects of a finite viscosity were 
taken into account. Nonetheless, the initial excitations 
are all in phase, as the interface starts with a flat profile, 
so ( [5.13 ) is a correct description of the short-time non- 
equilibrium growth of the interface width. 

From ( |5.10D and ( |5.13| ) one finds that W 2 scales accord- 
ing to 



W 2 (L,t) = L 3 g(t/L^ 2 ) , 
where for large times 

while the short-time limit is given by 
k B T 



(5.14) 



(5.15) 
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(t/L^ 2 ) 



ieir 4 



E—j sii 




(5.16) 



One can obtain an analytical form for the above expres- 
sion by appr oxima ting the summation with an integra- 
tion, so that ( 5.13 ) becomes 



W 2 {t) 



W 6/s 
57r(2p) 3 / 5 e 2 / 5 



dx 



sin 2 (x) 

.T n /5 



(5.17) 



where we set x — uj (k)t. The integral evaluates to |54 



dx 



sin 2 (a;) 5 2 2 1 / 5 



,11/5 



cos(2tt/5) T(4/5) = 1.7219. 



(5.18) 



and the roughening of the interface thus scales according 
to 



W 2 {t) 



2 3 / 5 5cos(27r/5)r(4/5) k B T 



3tt 



3/5 ,2/5 



P 



t 6 ? 5 . (5.19) 



We expect (5.19) to describe the growth of the interface 
until the crossover time in which equilibrium is attained. 
The crossover time t c may be approximated by the time 
necessary for the longest wavelength excitation to reach 
its maximum amp litude, that is t c = T /A. If the disper- 
sion relation (4^) is used to estimate the period T, we 
conclude that 



tc 



n 

2 w 



1 . pL* 



(5.20) 



We shall compare the lattice-Boltzmann numerical sim- 
ulations with the theoretical results of this section in the 
following. 



VI. THERMAL FLUCTUATIONS OF ELASTIC 
INTERFACES: SIMULATIONS 

In this section we present the results of the simulations 
performed using the lattice-Boltzmann method with a 9- 
velocity square latticejthat is, an FCHC lattice projected 
in 2-dimensions) . The computational box of size 

LxL is filled with two immiscible fluids pre-thermalized 
to a common temperature and separated by a thin inter- 
face. The two fluids have the same densities and viscosi- 
ties. We study the dynamical effects of thermal fluctu- 
ations at equilibrium and out of equilibrium, wh en th e 
interface grows up to the steady state given by ( 5.10 ). 
The interfaces we consider are driven by the coupling 
with the surrounding fluids and by purely elastic forces, 
tuned by the bending modulus e, as defined in Section y. 



A. Steady state 

The evolution of an interface of size L = 64 was mon- 
itored for 2 20 time steps. The resulting log-log plot of 
the frequency power-spectrum is given in Figure ^. We 
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FIG. 5. Frequency power-spectrum for a fluctuating inter- 
face of size L = 64, whose time evolution was monitored 
for 2 20 time steps. Experimental data for wavenumbers 
k — 2n/L,4n/L and 6n/L axe reported as filled circles, 
squares and diamonds respectiv ely. The solid lines represent 
the theoretical predictions from (5.2), showing that the agree- 



ment spans about three orders of magnitude. 



distinguish three regions: 

• A plateau corresponding to the low-frequency limit 

of id), 



lim \hk{u)\ /&= — 5 o— 

«-x>' 7T 2 (ek 4 ) 2 L 



(6.1) 



that shows explicitly a dependence on the bending stiff- 
ness e M . Low- frequency excitations are thus controlled 



8 



by elastic forces and viscous effects. In figure this region 
spans a frequency interval of about two orders of magni- 
tude, up to where u> approaches the pole of the response 
function T. 

• A transition region in proximity of the b end ing wave 
freq uen cy given by the dispersion relation (4.4), a pole 
for ( |5.2| ). Notice that a peak at such frequency manifests 
itself at higher wavenumbers. This behavior is opposite 
to the one observed for interfaces driven by surface ten- 
sion f||§. 

• The high-frequency region, expected from (5.2) to 



represent a u~ 7 / 2 decay, independent of either the bend- 
ing rigidity or the surface tension of the interface, accord- 
ing to 



lim \h k (u)\ 2 /e = 



8n 2 Lp U 



7/2 



(6.2) 



The relative amplitude of this segment of the power spec- 
trum is thus determined by viscosity effects alone. 

The experimental results reproduced in Figure || agree 
with the theoretical predictions, here drawn as solid lines, 
over a range of about three orders of magnitude, includ- 
ing part of the high-frequency region. The discrepancy 
observed at higher frequencies is probably due to stand- 
ing sound waves that create a distortion in the inter- 
face profile Since the theoretical prediction (5.2) 
has been confirmed by our numerical results, the subse- 
quent conclusions, equations (5.8) and (5.IC), describe 
the average wavenumber power spectrum and the func- 
tional dependence of the interface width on t he linear 
dimension of the box. In particular, from (5.10) we esti- 
mate that W/L « 1/16 for the case just discussed. 



B. Effective surface tension 

More simulations were performed to verify the agree- 
ment with the theoretically predictable wavenumber 
power spectrum. Different temperatures, lattice sizes and 
bending rigidities were prescribed to the system. In Fig- 
ure ^ we report a study of equilibrium power spectra for 
a lattice of size L = 60 and bending modulus e ranging 
from 0.0032 to 0.512. With this set of parameters the 
ratio W/L spans an interval from 1/6 to 1/60. We notice 
a very good agreement with ( |5.8| ) when W/L > 1/20, but 
for smaller ratios the relation describing the steady state 
power spectrum is approximated by 



(a cS k 2 + efc 4 ) L 



(6.3) 



showing the presence of an effective s urfa ce tension a e s. 
By considering the small-fc limit of (1x3) when W/L < 
1/20, we measured cr e ff for different points in the 3- 
dimensional e — L — T space. The following empirical 
scaling law resulted from our simulations: 



creff = ^V[(W/L) 2 ] = i V(TL/e) , 
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FIG. 6. Power-spectra of fluctuating membranes of size 
L = 60 lattice units and different bending stiffnesses. The 
interface- width/length ratio, W/L, ranges from 1/60 to 1/15. 
When W/L > 1/20 the purely elastic functional dependence 
of the power spectrum (with slope of —4 in a log-log plot) is 
recovered. 



where I? is a scaling function. Figure shows a log-log 
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FIG. 7. Effective surface tension as results from the power 
spectra of the simulations of fluctuating membranes. Com- 
putational boxes with L = 40 and 60, different temperatures 
and bending rigidities were considered. We report in figure a 
log-log plot of the data for a e a versus ksTL/e = W 2 /L 2 . 



plot of the measured values of er c ff versus (W/L) 2 . <r c s 
decays exponentially for increasing values of W/L; also, 
when W/L > 1/20 it practically vanishes and the purely 
elastic dynamics takes place. 

we 



L 2 



(6.4) 



In order to discuss the origin of equation (6.4), 
investigate here the scaling properties of the ratio 7/7, 
the crucial factor in the perturbation ( 3.16] ) of the Boltz- 
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mann equation. For a boundary of the form h(x), 
the derivative with respect to the arc length s can be 
expressed in terms of x as 



d 
ds 



l + h'(x) 



d_ 

2 dx ' 



(6.5) 



where the prime stands for dj dx. From the expression 
for the curvature, 



h"(x) 



(l + h' 2 ) 



3/2 



(6.6) 



one obtains 

7 ftW 



10- 



h'h" 



/2^7 1,2 



7 h"{l + h' 2 ) (1 + h' 2 ) 



t (l-5h u )h 
1 {l + h' 2 f 



(6.7) 



In the case of sinusoidal bending waves and fluctuating 
interfaces with horizontal periodicity as a boundary con- 
dition, the profile h(x) satisfies pa] 



W 



T n {2ttx/L) , 



(6.8) 



where h^ n \x) is the n-th derivative of h( x) a nd T n i s 
periodic in x with period L. By replacing (pTSJ) in (3.7), 
we find the following scaling relation for the dynamical 
factor, 



7 

7 



1 

L~ 2 



' Gth 



(W/Lf,x/L 



(6.9) 



where Gth is the theoretical scaling function. 

After initializing the system with different sinusoidal 
waves, we measured the ratio 7/7 for L — 64 and 256 
lattice units and amplitudes ranging from 1/64 to 1/8 
of the wavelengths. In Figure || we report the results of 
such measurements as a function of x. From the figure 
the experimental behavior of 7/7 presents the scaling 



7 

^ y ex 



(W/L) 2 ,x/L 



V 



(W/L) 2 ]}, 

(6.10) 



where Q ex approaches Gth for large computational boxes 
(high resolution): 



lim G e 



(W/L) 2 ,x/L 



Gth 



{W/L) 2 ,x/L . (6.11) 



The most notable feature in ( |3.1C| ) is the presence of the 
additional constant term T>[(W / L) 2 ] / L 2 . It is responsible 
for the introduction in our model of an effective surface 



tension, as it causes S in (3.16) to be a non- vanishing 
number. 
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FIG. 8. Measured values of the ratio 7/7 are plotted as 
filled circles for static sinusoidal interfaces with wavelengths 
L of 64 and 256 lattice units and widths W ranging from 1/64 
to 1/8 of the wavelengths. Solid lines represent the anal ytica l 
results for the continuum limit obtained from equation (6.7). 
When W/L < 1/32, the measured values appear to be shifted 
with respect to the theoretical ones. This introduces an effec- 
tive surface tension for almost flat surfaces that scales as a 
function T> of W and L: a cB = V[(W/L) 2 ] / L 2 . 



C. Non-equilibrium roughening 



As shown in Section VB, at large times the average 
i nterf ace width reaches its equilibrium value given by 
( |5.10 ). Before reaching the stationary state, thoug h, th e 
interface growth is described by the power-law (5.19). 
In this section we present the results of the simulation 
of non-equilibrium growing interfaces and compare them 
with the theoretical predictions. 

After initializing the system as a flat boundary separat- 
ing two immiscible pre-thcrmalizcd fluids, we monitored 
the growth of the interface for a period of time approx- 
imately equal to twice the crossover time t c estimated 
in ( 5.20p . We repeated the simulation a hundred times 
and averaged the results, in order to obtain the mean 



quantities that appear in equations (5.10) and (5.19) 



Such an experiment is reported in Figure]^, where the 
solid line represents the expected behavior of the inter- 
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FIG. 9. Log-log plot of the time evolution (circles) of the 
mean width of an interface with size L = 64, starting from 
flatness at time t = 0, and monitored fo r an interval of time 
equal to about twice the crossover time ( 5.2C ). The solid line 
represents the theoretical scaling law W oc t / 3 predicted for 
the effective surface tension <r ff. The sudden change in slope 
around t ~ 10 3 corresponds to the onset of the purely elastic 
dynamics. 



face growth, taking into account the presence of the effec- 
tive surface tension. From the figure, the early-time evo- 
lution of the interface growth appears to be described 
by the scaling law W oc i 1 / 3 |52|. Again, we conclude 
that the membrane is indeed driven by a c s for very 
small interface- width/length ratios, that is for almost fiat 
boundaries. In order to by-pass the effect of the surface 
tension and analyze the purely elastic non-equilibrium 
dynamics, the membrane profile is initialized in such a 
way that W/L > 1/20. Under these conditions, the dis- 
turbance due to cr e ff is negligible (as shown in Section 
VIA). The Fourier components of the interface profile 
are prescribed as 



hk(t ) = y ~~7T4~ sm[LU (k)t \ , 

where t Q — t c /A is the time at which the numerical 
simulation starts; furthermore, the relative intensities 
of hk{t ) have been chosen so that they agree with the 
evolution (5.12). We studied the non-equilibrium rough- 
ening of such interfaces for an interval of time starting 
at t c /A and terminating at approximately 3t c , after the 
equilibrium width is attained. In Figure frCl(a) we report 
the interface widths for system sizes L = 32,48,64 and 
96, rescaled as W/L?/ 2 so as to match the prefactor 



^/fcsT/720e from the anticipated steady state (5.10), 
drawn in figure as a solid line. The horizontal coordi- 
nate in Figure |lQ(a) represents the time t in lattice units, 
rescaled as t/Lr> 2 so that all the crossover times fall at 
y/p/Mir 3 e, according to (5.20). In order to evaluate the 
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FIG. 10. In (a) we plot the rescaled interface width W/L 3/2 
in lattice units, as a function of the time, scaled to i/L 5 / 2 
so as to make all the cross over times match the prefactor 
yp/647r 3 e, predicted from (5.2C). The computational box 
sizes are L = 32, 48, 64 and 96. The horizontal line represents 
the steady state width as in (5.10). In (b) a log- log plot of 
the same data is given. The dashed line marks the scaled 



crossover time, while the obliq ue lin e represents the theoret- 
ical prediction from equation ( |5.19| ), with an angular coeffi- 
cient of 3/5. The symbols represent different values of the 
system size: squares, diamonds, filled circles and six-pointed 
stars correspond to L — 32, 48, 64 and 96 respectively. 



prefactors above, we used the relations fl3.15| ), ( p.l7[ ) and 
(4.5), where the parameters for this experiment were set 
to Xb = —1.5, E = 0.0002, p = 0.5 and the variance A to 
10 -4 . A log-log plot of the same data is sh owed in Fig- 
ure |o|(b). Also, the predicted scaling law ( 5.19| ) is rep- 
resented as a solid line with angular coefficient 3/5. We 
notice a good agreement between the numerical results 
and the theoretical predictions, although we were able to 
investigate the non-equilibrium behavior in a time range 
of just about one-half order of magnitude. 
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VII. CONCLUSIONS 

Fluctuating elastic interfaces in fluids have been dis- 
cussed from different approaches. By means of the 
fluctuation-dissipation theorem, we presented a theoret- 
ical derivation of the non-equilibrium interface roughen- 
ing, described by the scaling law W 2 = L 3 g(t/L 5 / 2 ) . 
The above expression predicts that the interface, initially 
flat, grows as W oc i 3 / 5 , and at long times it reaches 
the equilibrium width W = ^k B T/720eL 3 / 2 . Also, we 
showed that the correlation function of the fluctuations 
decays exponentially in the long-time limit, in agreement 
with the behavior of the correlated forces on the inter- 
face, due to the excitations in the bulk ]52| . 

Although the theoretical discussion was necessary for 
the sake of clarity, this work instead focused on the 
development of a numerical method simulating the phe- 
nomenology of fluctuating membranes. We created a 
lattice-Boltzmann method. Starting from the descrip- 
tion of the membrane dynamics by means of the Landau- 
Helfrich free energy |I],[mJ], in which the interface geome- 
try appears explicitly in the form of the curvature 7 and 
its derivatives, we translated the macroscopic equation of 
motion governing the interface evolution in a perturba- 
tion of the single relaxation-time lattice-Boltzmann equa- 
tion. Such a perturbation depends on 7, so that a crucial 
point of our model is the localization of the interface and 
the measurement of its geometric properties. As outlined 
in Appendix [B|, we adopted an explicit characterization 
of the boundary followed by a polynomial mapping, from 
which we extracted the information about the local cur- 
vature. 

We tested our method to reproduce the hydrodynamic 
equilibrium of circular bubbles and the dynamical cou- 
pling of the interface with the surrounding fluids in the 
bending wave dispersion relation. Thermal fluctuations 
were introduced in the model by adding a random compo- 
nent to the fluid stress tensor |4^] . The lattice-Boltzmann 
equation was thus generalized to include a stochastic 
term 35 1, whose fluctuations are uncorrelated in space 
and time. The equilibrium frequency power spectrum 
of fluctuating elastic interfaces, predicted by the theory 
and related to the correlation function of the excitations, 
was confirmed in our numerical simulations. Also, we 
simulated the non-equilibrium roughening of membranes, 
monitoring the time evolution of the interface width. Its 
growth rate resulted in agreement with the theoretical 
predictions mentioned above, although we noticed a dis- 
turbance due to numerical errors at almost flat surfaces, 
generating an unwanted effective surface tension. In 
order to improve the effectiveness of the interface detec- 
tion process, other, more accurate, techniques for track- 
ing the interface could be considered. For example, using 
markers in the context of an explicit discretization of the 
interface J57],^8| could further improve the performance 
of this numerical method. 

One of the advantages of our numerical method relies 



on the fact that 7 and d 2 j/ds 2 , the second derivative in 
the arc length s, are not approximated by their linearized 
expressions d 2 h/dx 2 and d 4 h/dx 4 respectively (here h(x) 
is the boundary profile). We could thus use our model 
to study and simulate systems in the nonlinear regime, 
that is when the interface is far from being flat and the 
curvature is a nonlinear function of the interface pro- 
file. Another feature of this method is that it has been 
designed from the outset to describe thin interfaces sep- 
arating two (or more) fluids, a situation that is diffi- 
cult to describe in terms of a slowly-varying continuous 
field (order parameter). In future work, we would like 
to employ this model in the study of fluctuating mem- 
branes separated in many distinct components, colliding 
with each other and immersed in fluids with prescribed 
or complex fields. Also interesting would be the inclusion 
of a surfactant species, as in the microemulsion model of 
Ref. [El. 



ACKNOWLEDGMENTS 

We thank Tom Chou and Michael Brenner for valu- 
able discussions. We are also indebted to Tony Dins- 
more, Peter Sheridan Dodds and Joshua Weitz for useful 
comments. 



APPENDIX A: FLUID-INTERFACE COUPLING 

In the following we review the dynamical coupling 
between the membrane and the surrounding fluid, gen- 
eralizing the results to rt-dimensional interfaces. The 2D 
case with a saddle-splay term in the free energy (2.1) is 
also presented. We shall denote with r = (x±, . . . , x n +i) 
the locus of the n-dimensional hypersurface imbedded 
in the (n + l)-dimensional space. The hypersurface is 
parameterized by the n coordinates u\, . . . , u n , collec- 
tively referred to as u, so that r = r(u). n(u) represents 
the unit normal to the surface oriented outwards and 
e; = dr/dui is the tangent vector corresponding to the 



coordinate u t . We shall assumes that the hypersurface is 
well-behaved and we are able to choose a parameteriza- 
tion u such that 

{e h e j )=e i -e j =8 ij . (Al) 

The Weingarten operator, d efin ed as L(ei) = dn/dui, is 
self-adjoint with respect to ( |Al[ ), and the e,'s are chosen 
to be eigenvectors of L 



L(et) 



dn 

dui 



7> e * 



(A2) 



where 7, is the principal curvature corresponding to e;. 
In other words, our parameterization u runs along the 
principal directions of curvature. We shall apply Hamil- 
ton's principle fl59| to the line integral of the Lagrangian 
describing the interface dynamics 
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lid = 



lr % -%H 2 -a)dSdt. 



(A3) 



where p is the density of the hypersurface, e is the bend- 
ing rigidity, H = 7^ is the mean curvature and a the 
surface tension. dS — du\ . . . du n is the volum e ele ment 
for the given parameterization. By means of (Al) and 
(1A2) one can verify that 



d 2 r 
du 2 



-an, 



(A4) 



so that (Kq) is recast as 



P -r 2 



lid — 

where 

dr 

A variation of ([ 

Slid = 



{V 2 u r) 2 -a]J,\d Ui r\ \ dSdt (A5) 



V 2 r = 

v u 

3) reads 



^ duV 



\du,r\ 



pr ■ -^Sr - eV 2 r ■ V 2 u 5r 



t-r 1 oui oui 1 



(A6) 



Upon integration by parts, the boundary terms vanish 
as we study closed surfaces or impose periodic bound- 
ary conditions and all of the derivatives involved are 
supposed to be smooth functions of the coordinates u. 
Therefore, by applying Hamilton's principle, including a 
term for an external force F ext 



■ SrdSdt 



(A7) 



5I = Jj[-pr-e{Viyr + aVlr 
+ J F ext ■ drdt = , 
one retrieves the following equation of motion 

F ext = -aV 2 u r + e(Vl) 2 r + pr . (A8) 

If F ext is the pressure, it is directed along the normal to 
the surface, that is F ext = Fn. By noticing that (A4) 



and (A2) lead to 



■r = (V 2 u )(-Hn) = 



OH 
< in , 



HP 1 

OUi 



(A9) 



and by projecting r along the normal and tangential com- 
ponents, (AS) can be separated as 



F — pr ■ n = aH - eV 2 H + eH £\ 7l 2 
pr ■ ei/e = 2 7i dH/dui + Hdji/dui 



(A10) 
l...n (All) 



where we have dropped the subscript u in the Laplacian 
as it is invariant for any isometric reparameterization. 



The main result of this Appendix is that the dynami- 
cal coupling between the interface and the flui d in the 
lattice-Boltzmann method is provided by ( AlCj ). Since 
the macroscopical motion of the membrane is extremely 
slow in terms of lattice units, of O(10 3 ) time steps, 



the second term in the left-hand side of ( |A10| ) can be 
neglected for our purposes and we shall regard 



F 



aH-e\7 2 H + ei/V 7 2 



(A12) 



as a quasi-stationary equilibrium between the interface 
configuration and the surrounding fluid pressure. Simi- 
larly it can be showed that the dynamical coupling for a 
2D interface, including as in (2.1) a saddle-splay correc- 
tion to (A3), is given by 



F = aH — e\7 2 H + eH(H 2 - 2K) 
I (HK - d 2 11 /du 2 2 - d 2 l2 /du 2 1 ) . 



(A13) 



where K = 7172 is the Gaussian curvature. ( A13| ) should 
be used in simulating 2D membranes with this lattice- 
Boltzmann method, when changes in the interface genus 
are taken into account. 



APPENDIX B: EVALUATION OF THE 
CURVATURE 

In this appendix we include a description of the algo- 
rithm used in this work to measure the geometric prop- 
erties of the interface. The procedure is schematically 
divided into three steps: localization of the interface, 
polynomial approximation, and evaluation of the local 
curvature. 



As mentioned in Section III, the populations corre- 
sponding to two immiscible fluids are distinguished from 
each other by splitting the occupation numbers n,(x, t) at 
a given lattice site in a "red" part, rj(x, t), and a "b lue" 
part, &i(x, t), of the distribution function, as in (3.8). By 
defining the color field as 



*(x,t) 



E 



ri(x,t) - &i(x,i)] 



(Bl) 



one can visualize the binary fluid as a 2D surface where its 
lowest regions correspond to the physical presence of the 
blue fluid and the flat "highlands" to the areas occupied 
by the red fluid. The interface, in this picture, is repre- 
sented by the steep slopes. One can define the boundary 
B between the two fluids as the set of lattice sites whose 
color field absolute value is smaller than a fraction a 
(< 1) of the fluid density p, that is 



£={x: |$(x,t)| <ap} 



(B2) 



In our simulations, due to the small thickness of the inter- 
face, a might be chosen to be a value between 0.5 and 
0.9 with no substantial changes in B. Another param- 
eter used in the description of the interface is the color 
gradient f , defined as 
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f(x,t) 



[r i (x,i)-& i (x,t)]. 



(B3) 



From its definition f appears to be, in the lattice approx- 
imation, perpendicular to the interface. The vector 
n = f/|f| will be used here as the unit normal to the 
interface. 

After localizing the interface, in order to measure its 
geometrical properties at a certain lattice site x Q in B, 
we replace a neighborhood of x Q with the polynomial 
resulting from a least-square fit, as described below. The 
neighborhood of a given boundary point (in the present 
case x ) is constituted by contiguous lattice sites in B 
and is chosen so that the following two conditions are 
met: 

• The normals at the points furthest from x D make an 
angle of at most 7r/4 with the normal at x D . 

• All the boundary points in the neighborhood lie 
within a distance of 30 lattice units from x D (about 10 
times the interface width) or half the linear size of the 
computational box, whichever is smaller. 

In Figure [n] we present the two-step process. The 
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FIG. 11. Filled circles represent an interface neighborhood 
of a given lattice site x D on the boundary (marked by a cross). 
An arrow indicates the direction of the normal n to the inter- 
face at the corresponding lattice site. The whole set of points 
in (a) is rotated clockwise by the angle ■& — arccos(n ■ x) so 
as to align n with the t/-axis (b). A least-square fit of the 
interface can thus be realized by the single- valued polynomial 
y = J2t=o aiX% /*' ' nere drawn as a dashed line. 



filled circles represent the lattice sites belonging to the 
boundary in a neighborhood of the point x G (marked by 
a cross) ; arrows are magnified representations of the nor- 
mals at the corresponding lattice sites. As shown in Fig- 
ure |ll|(b), calculations are simplified by rotating clock- 
wise the whole set of points by the angle z? = arccos(n- i), 
where i is the unit vector directed along the positive x- 
axis. In this way, we can fit the subset of B by a single- 
valued polynomial in i: y = YliLo o-%x % /i\ , where x = 
corresponds to the horizontal coordinate of x Q . It is con- 
venient to choose N > 4, because derivatives up to the 
fourth order are involved in the evaluation of 7, the sec- 
ond derivative in the arc length of the curvature. In 
practice, we let N = 5, as the accuracy of our method 
does not improve for larger values of N. 



Once the coefficients a, are estimated, we c alcula te 7 
and 7/7 by replacing h{x) respectively in (6.6) and (6.7) 
with y = ^2 ciiX 1 /i\ and set ting x to van ish. In other 
words, a n substitutes hy 1 ' in (3.6) and (6.7). The dynam- 
ical term 7 2 (x, t) — 7(x, i)/7(x, t), appearing in the per- 
turbation ( 3.16| ), is thus evaluated by repeating the above 
procedure for each lattice site x in B. 
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